#ifndef _MATRIX_H
#define _MATRIX_H
#include<bits/stdc++.h>
using namespace std;

template <class T> class matrix;
template <class T> class rowvec;
template <class T> class colvec;

template <class T> T operator * (const rowvec<T>&,const colvec<T>&);

template <class T> class matrix{
private:
	int n, m;
	T *a;
public:
	matrix(){
		a = NULL, n = m = 0;
	}
	matrix(int _n,int _m){
		if(_n <= 0 || _m <= 0){
			cout << "Illegal Creating!" << endl;
			exit(-1);
		}
		n = _n, m = _m;
		a = new T[_n * _m];
		memset(a, 0, sizeof(T) * _n * _m);
	}
	const T& at(const int& i)const{return a[i];}
	T& at(const int& i){return a[i];}
	const T* operator [] (const int& i)const{
		return &a[i * m];
	}
	T* operator [] (const int& i){
		return &a[i * m];
	}
	pair <int,int> size()const{
		return make_pair(n, m);
	}
	bool issq(){
		return n == m;
	}
	void read(){
		for(int i = 0; i < n * m; ++i) cin >> a[i];
	}
	void write()const{
		for(int i = 0; i < n; ++i, cout << endl)
			for(int j = 0; j < m; ++j)
				cout << a[i * m + j] << " ";
	}
	matrix(const matrix& p){
		n = p.n, m = p.m;
		if(n && m) a = new T[n * m];
		for(int i = 0; i < n; ++i)
			for(int j = 0; j < m; ++j)
				a[i * m + j] = p[i][j];
	}
	const matrix& operator = (const matrix& p){
		n = p.n, m = p.m;
		if(a != NULL) delete a;
		if(n && m) a = new T[n * m];
		for(int i = 0; i < n; ++i)
			for(int j = 0; j < m; ++j)
				a[i * m + j] = p[i][j];
		return *this;
	}
	~matrix(){
		if(n && m && a != NULL) delete a;
	}
	
	friend matrix<T> operator * (const matrix<T>& a, const matrix<T>& b){
		if(a.m != b.n){
			cout << "Error: Illegal Multiply!"<<endl;
			exit(-2);
		}
		int n = a.n, p = a.m, m = b.m;
		matrix <T> c = matrix <T> (n,m);
		for(int k = 0; k < p; ++k)
			for(int i = 0; i < n; ++i)
				for(int j = 0; j < m; ++j)
					c[i][j] += a[i][k] * b[k][j];
		return c;
	}
};

template <class T> class rowvec : public matrix<T>{
public:
	rowvec():matrix<T>(){}
    rowvec(int _n):matrix<T>(1, _n){}
	const T operator [] (const int& i)const{
		return matrix<T>::at(i);
	}
	T& operator [] (const int& i){
		return matrix<T>::at(i);
	}
	int size()const{
		return matrix<T>::size().second;
	}
	friend T operator * (const rowvec<T>& a,const rowvec<T>& b){
		if(a.size() != b.size()){
			cout << "Illegal Inner Multiply!"<<endl;
			exit(-3);
		}
		T res = 0;
		int n = a.size();
		for(int i = 0; i < n; ++i) res += a[i] * b[i];
		return res;
	}
	friend T operator * <> (const rowvec<T>&, const colvec<T>&);
};

template <class T> class colvec : public matrix<T>{
public:
	colvec():matrix<T>(){}
	colvec(int _n):matrix<T>(_n,1){}
	T& operator [] (const int& i){
		return matrix <T> ::at(i);
	}
	T operator [] (const int& i)const{
		return matrix <T> ::at(i);
	}
	int size()const{
		return matrix <T> ::size().first;
	}
	friend T operator * (const colvec<T>& a, const colvec<T>& b){
		if(a.size() != b.size()){
			cout << "Illegal Inner Multiply!"<<endl;
			exit(-4);
		}
		T res = 0;
		int n = a.size();
		for(int i = 0; i < n; ++i) res += a[i] * b[i];
		return res;
	}
	friend T operator * <>(const rowvec<T>&, const colvec<T>&);
};

template <class T> T operator * (const rowvec<T>& a,const colvec<T>& b){
	if(a.size() != b.size()){
		cout << "Illegal Multiply!"<<endl;
		exit(-5);
	}
	T res=0;
	int n = a.size();
	for(int i = 0; i < n; ++i) res += a[i] * b[i];
	return res;
}

#endif